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A gas-kinetic scheme for the continuum regime is apphed to the simulation of turbu- 
lent compressible flow, by replacing the molecular relaxation time with a turbulent 
relaxation time in the BGK model. The turbulence dynamics is modelled on the basis 
of a standard, linear two-equation turbulence model. The hydrodynamic limit of the 
resulting turbulence model is linear in smooth flow and non-linear in the presence 
of stronger flow gradients. The non-linear correction terms in the numerical flux are 
weighed as a function of "rarefaction" - referred to turbulence dynamics and not to 
molecular dynamics, i.e. measured by the ratio of turbulence to mean flow scales 
of motion. Even though no assumptions on the nature of the turbulence have been 
made and a linear two-equation turbulence model is used, the turbulence gas-kinetic 
scheme seems able to correct the turbulent stress tensor in an effective way; on the 
basis of a number of turbulence modelling benchmark flow cases, characterized by 
strong shock - boundary layer interactions, the turbulent gas-kinetic scheme provides 
a level of accuracy comparable to that of more sophisticated turbulence models. 
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I. INTRODUCTION 



Gas-kinetic schemes (GKS) simulate the mechanics of fluids by modelling the numerical 
fluxes between computational cells on the basis of the gas kinetic theory - often relying 
on the BGK modeP. GKS have been developed for the continuum regime with the aim 
of achieving a better accuracy and stability than traditional numerical schemes based on 
the Navier-Stokes equations (NSS). GKS are consistent with traditional schemes as the 
Euler and the Navier-Stokes equations can be derived both from the Boltzmann equation 
and from the BGK model, for instance by means of the Chapman-Enskog expansiorpEl. A 
number of GKS were developed over the past twenty year^3HHl_ Among them, the scheme 
developed by Xu in 2001^ has achieved a signiflcant level of validation, having been tested 
in a number of laminar flow cases in the continuum regime, ranging from low-Reynolds 
subsonic to hypersonic^^'^^. The scheme has also been tested in the transitional regime. 

The results provided by this GKS show that it disposes of mechanisms able to recover the 
gas-kinetics physics and is potentially able to provide physically more consistent solutions 
than NSS. In practice: it provides (i) physically consistent solutions in all tested flow regimes, 
(ii) solutions very close to those provided by NSS in smooth flow conditions, (iii) accurate 
solutions also inside shock-layers, which differ from NSS. The higher accuracy of the GKS is 
attributed to the underlying gas-kinetics which dictates a time-space coupled gas-evolution, 
whereas the numerical fluxes in NSS are time-independent. Moreover, advective and viscous 
fluxes are calculated simultaneously, whereas in NSS the advective fluxes are usually provided 
by the solution of the Riemann problem, which does not consider the viscosity of the fluid 
or the effects of the particle collisions. In laminar flow, the effects of collisions on the 
advective fluxes are only signiflcant in shock layers, where the shock thickness is comparable 
to the particles mean free path. In turbulent flow, the rarefaction measured on the basis of 
turbulence dynamics, is much higher: a turbulent Knudsen number calculated on the basis 
of eddies mean free path, can locally reach values well in excess of the 0.001 threshold. 

This study proposes an analysis of the behaviour of a GKS - the one developed by XiP - 
applied to the simulation of turbulent flows. The development of a turbulent GKS is straight 
forward, provided a consistent eddy viscosity fleld is available: the molecular relaxation time 
in the BGK model is replaced by a relaxation time consistent with the dynamics of turbulent 
structures. However, the implications of this replacement are more profound than it may 
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seem. According to Chenl^^'^ this replacement implies the modelling of a distribution of 
eddies instead of a distribution of molecules. It follows that turbulence is modelled by 
a "mixing time" instead of a "mixing length", and that the hydrodynamic limit of the 
"turbulent" BGK model exhibits a non-linear turbulent stress tensor, whose mathematical 
structure is identical to that of non-linear turbulence models. 

A number of flow cases have been simulated with the GKS and a solver based on a 
standard NSS. Turbulence has been modelled in both solvers with a standard, linear k- 
uj model. This class of turbulence models, when used in NSS, fails to accurately position 
separations points and tends to underestimate the size of separation areas. The selected flow 
cases include strong shock - boundary layer interaction; in order to test the GKS behaviour, 
they have been selected among those where standard NSS fail to provide accurate solutions 
with linear turbulence models. 

This paper is structured as follows: the GKS is presented in section [Tl| the turbulent 



fluxes generated from a turbulent relaxation time are analysed in section III , the analysis of 



the turbulent GKS is presented in section IV, the numerical implementation is described in 
section |V| the numerical experiments are presented in section |VI[ conclusions are exposed 
in section IVIII 



II. GAS-KINETIC SCHEME 

Following Xu's papei'^ the scheme assumes a finite-volume discretisation with discon- 
tinuous reconstruction of the conservative variables. The macroscopic variables are w = 
[p pui pu2 pus pE]'^, where p, U, and E indicate density, velocity and total energy of a gas. 
The BGK model^ is used instead of the Navier-Stokes equations: 

df ^ df jg-f) 

where the summation convention holds, / is the gas distribution function, g is the equilibrium 
state, a Maxwellian distribution, approached by /, and r is the relaxation time, taken to be 
the time between collisions, which is related to the molecular viscosity and heath coefficients 
of the gas r = p/p. The equilibrium distribution is: 



9 = P 



K+2 



TC 



,-x{{u,-Uif+e) ^ (2) 

where A = m/ {2kT), m is the molecular mass, k is the Boltzmann constant, and T is 
temperature. The variable ^ indicates the N effective degrees of freedom of the gas molecules 
= (5 — 37) / (7 — 1) + 1, where 7 is the specific heat ratio. The macroscopic variables are 
obtained by taking moments of /: 



w = I ijfdE, (3) 

where in three dimensions is: 



1 ^ 

1 Ml U2 U3 - (Ui^ + + + 



(4) 



Conservation of mass, momentum and energy during particle collision is expressed by: 



{g-f)^dE = 0. (5) 
Eq. [T]can be readily solved analytically: 



1 f 

f{x, y, z, t, u, v,w,C) = - / g{x', y', z', t, u, v, w, Oe~^*~*'^/^ dt' +e'^^^ foix-ut, y-vt, z-wt), 

^ Jo 

(6) 

where /o is the initial gas distribution function, x' = x — u{t — t'), y' = y — v{t — t'), z' = 
z-w{t-t'). 

The numerical fluxes between two cells, assuming an interface normal to direction xi 
takes the form (per unit length): 

J'i = f f u.^PfdEdt, (7) 



Jo J 

where Ui and if) are evaluated at the interface. The distribution functions g and /o in Eq. 
[6] reflect the discontinuous nature of the assumed reconstruction and have a discontinuity 
xi = 0. The equilibrium distribution is expanded in all spatial directions and time up to 
the first order: 

^ _ go{l + d[x^ - At) , xi < 0, 
go (1 + d^Xi - At) , xi > 0, 
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where the superscripts / and r indicate the left and right states and go is a MaxweUian 
corresponding to a state Wq = [po po'^oi Po^02 Po'^os PoEq], which is an average between 
left and right states. The initial distribution is expanded in all spatial directions and time 
including non-equilibrium terms, the ones modelling the effects of collisions, i.e. the terms 
proportional to the relaxation time r: 

_ g\l + a\xi) - T [aim + A^) , xi < 0, 

Jo — \ (.yj 

5(''(l + a[xi) -r(a[Ui + A"), xi > 0, 
The choice of the terms used in the expansions in Eqs. [8] and [9] directly affects the 
accuracy of the scheme. The Taylor expansions in Eqs. [8] and [9] are multidimensional, 
as they include all coordinate directions. Some applications of GKS use the directional 
splitting^'^, where the coordinate x„ identifies the direction normal to 

the interface, based on the experimental evidence that tangential expansion terms provide 
negligible contributions to the fluxes. 

Each of the coefficients in Eqs. [8] and |9] is expanded in all degrees of freedom: 

ai = an + ai2U + a^gt; + + ai^{u^ + f ^ + + = [ai a^s an ^is]^ • (10) 

The expansion coefficients are determined by derivation in all coordinate directions from 
Eq-i 

' /■^.,dE=/'AE=<!'«"^'-=°' (11) 



dxj dxj I J dx 



■'I 



Mai, > 0, 



' (12) 

Mai, > 0, 



dxi dxi J J dx_ 

where the matrix M, which can be easily inverted'*, is 



M = j iP^^gdE. (13) 

The temporal expansion coefficients A and A are determined from Eq. [5]and by imposing 
that the corresponding flux is also zero: 

At 

'ig-f)^dE = 0. (14) 



The equilibrium state wo is also obtained from Eq. [5] 



wo = j i^g^h^dE + j ijg'^h^dE. (15) 

Once all expansion coefficients have been determined, the distribution / is obtained 
inserting Eqs. [8] and M into Eq. [6] 



/ = (1 - e-*/") go + (-r + re"*/" + t e~''^) {h' a} + K a/') mgo^it-r^ re"*/") Ag^ + 
+ e-*/" (/i'^' + /i"^" - (t + r) (uia^/i'^' + u,a\K g')) - Te'^l^ [A^h^g^ + A'ffg^) , (16) 

where H is the Heaviside function, W = H{U), = 1 — H{u) and all coefficients are 



intended as series expansions of the form introduced in Eq. [TO 

Two modifications are finally suggested in the original papei'^ (i) the BGK model implies 
a unity Prandtl number; the heath fiux must therefore be corrected for realistic gas / fiuids, 
(ii) it is suggested that the simulations are conducted with a modified relaxation time r, 
which includes an artificial dissipation term, proportional to the pressure jump across the 
interface, to the molecular viscosity. 

r^t + \f::J!iM. (17) 

p Ip''" + p I 

III. TURBULENT GAS-KINETIC SCHEMES 

The GKS presented in the previous section |TT] is modified in order to account for the effect 
of unresolved turbulence on the resolved mean fiow, following the approach often referred to 
as RANS in its implementation with NSS. The implementation of a turbulence model into a 
gas-kinetic scheme is straight-forward; a relaxation time is introduced, which considers both 
molecular and turbulent phenomena: 

r = (18) 

P 

In practical calculations, the relaxation time is corrected by adding an artificial dissipation 



term, as was done in the laminar case in Eq. 17 



r = ^ + K^At. (19) 



p \p^' + I 
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The turbulent relaxation time can also be expressed as a non-linear function of the macro- 
scopic variables as suggested by Chen (adapted from Chen's original formulatiorP^t) . Tur- 
bulent and total relaxation times can be expressed as: 



r = ^+ ''^\,, +c\f^At, (20) 

where t] = S/u is the ratio between unresolved and resolved turbulence time scales, where 
is a scalar representing local velocity gradient and the coefficient C is determined heuris- 
tically (with an upper limit in the order of one). Numerical fluxes are obtained from Eq. 

m 

By neglecting variations of eddy viscosity across the interface, the analytical solution 
shown in Eq. [6] is still valid. The distribution function / and the numerical fluxes are 



obtained just like in the laminar case, inserting Eq. [19] or Eq. |20]into Eq. [16 

Turbulent quantities and eddy viscosity are modelled with a standard linear k-oo modeP^ 
two-equation models, which is a well-known member of this class of models and is represented 
by the following equations: 

dpK ^ dpUjK ^ o* . , ^ ft , _* 



dt dxj K dxj \ dxj 



/if = 7 . (23) 

where P = Tijdui/dxj is turbulence production term: 
The turbulent stress tensor is given by: 

/ 2 duk ^ \ 2 , - _ ^, 

= /if \2Sij - 1 - ^pkdij, (24) 

and Sij = 1/2 {dui/dxj + duj/dxi) is the strain rate. 

A typical choice for the parameters is (3 = 3/40, (3* = 9/100, 7 = 5/9, 7* = 1, cr = 1/2, 
a* = 1/2. A number of different versions of the k-u model have been developed: this study 
is based on the 2006 version. In the turbulent GKS, the isotropic part of the turbulent stress 
tensor in Eq. [24] can be added explicitly to the fluxes computed with Eq. [7] 



IV. ANALYSIS OF THE RESULTING TURBULENCE MODEL 



The fluxes in direction Xi obtained from Eq. [7] are recast as a series expansion in r 
truncated to tlie second order: 

^ = ^(0)+ ^(1)^ + ^(2) (25) 

Tfie superscripts refer to the order of the term, for the sake of clarity the subscript i is 
dropped. The three expansion terms are given by: 

jr(o) = y ^ ^At + ^^t^A^ goUjdE (26) 
J^(^) = j tp[At [h^ aim + a{ui - A) + (1 - e) d + /i'^' - ^o) + ^ttDu^ Ujd^7) 
JF(2) = (1 - e) j ^ [{Ago - A^h^g^ - A'h'^g'') uj - 2Du, uj] dE (28) 

where e = e"^*/"^ and D = h''d/gQ + h''di'gQ — a\h'-g'- — a'^-h'^g'^. 

Eq. |25] is closely related to a Chapman-Enskog expansion, observing that the relaxation 
time T can be related to a non-dimensional quantity providing a measure of rarefaction in 
the "distribution of eddies", i.e. a turbulent Knudsen number Knt. It may be defined as the 
ratio between the eddies' mean free path, A(, and the mean flow lengthscale Im, Krit = Xt/lm- 
The eddies mean free path is estimated on the basis of the turbulent relaxation time Tt and 
the mean velocity of the flow U: 

Xt = 2TtU (29) 
whereas the mean flow lengthscale can be obtained from the gradient of density: 

max{dp/dxi) max (91n (p) /(9xj) ' ^ ^ 

leading to the following turbulent Knudsen number: 

Knt = 2^^, 

max {d In (p) / dxi) 

In order to provide some insight into the nature of the GKS fluxes obtained, that is, 
to characterize the hydrodynamic limit of the turbulent GKS scheme, Eq. [25] must be 
simplified. It can be represented as a sum of terms, each of which identifies a particular flux 
- either advective {Fa) or viscous (Fy) - pertaining to the equilibrium state (no symbol) 
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or expressed as a weighted average between left and right state (indicated with the tilde ~ 
symbol). Furthermore, the symbol indicates time-accurate fluxes, which include the time 
coefficient A or A. Introducing the following time-independent fluxes: 



-pta 



-pta 



- J ip [h'' d/ui + h'' di'ui + A) goUidE, 

ip {a\uih!'g^ + a^Uih^g^) uidE, 
J ip [h'- aiuig^ + h"" ai'uig'' + A) uidE, 



(32) 
(33) 
(34) 
(35) 
(36) 
(37) 



Eqs. p6, 27 and 28 can be further simplified. The fluxes can be reexpressed as follows: 



1 . 

2^' dt ' 



U 7a- T A 



+ ^ ^CFL Ax, 

J-(2) = (l-e)(j-^-J-*? + J-v.-Jv 



(38) 



(39) 



CFL 



Axi 



Axi 



(40) 



where the relation At = CFLAxi/U has been used, CFL is the Courant-Friedrichs-Levi 
number and Asi the size of the computational cell in direction xi. d7A/dt = J ipAgoUidE, 
is the time derivative 7a- 



The finite differences in Eqs. |38} [39] and |40] have no physical meaning, but they pro- 
vide additional insight as each of these differences can be related to the derivative of the 
corresponding flux in direction Xi. For instance. 



— d7v . 

J-y — J-y = ay— Axi 

OX I 



(41) 



where ay is a coefficient dependent on the local numerical reconstruction. Similarly, coeffi- 



cients aA and ay are introduced. Eqs. 39 and 40 can then be re-expressed as: 



^<...^.(^.,,„,^^,.,,_^„,^). (42) 

JFP' = Ai(l-6)-|7(a,~+av)?P^. (43) 

The fluxes' derivatives introduce non-linear terms. The structure of these non-linear 
terms can be explained by considering the simple case of smooth flow, where right and left 
interpolation are identical; in this case the viscous fluxes generated by the GKS can be 
expressed as: 



Tv = - \ ipttiUi goUidE. (44) 



Recalling Eq. 11 the coefficients Oj can be expressed as = M dw/dxi. Eq. 
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can 



therefore be rewritten in the following form: 

Ty = NaM-'^^, (45) 

where the matrix Nn is given by the integral Nn = — J ip goUiUidS. The derivative of the 
viscous fluxes can then be reexpressed: 

^ i AT «/f— i ^1 AT I AT «/r— i 



dxi dxi 



( N„ l^x.) = A M-) + N„ ^. (46) 

\ OXi J OXi ^ OXi oxi 



Eq. 46 shows that the high-order terms in GKS fluxes are non-linear and depending on 
the numerical reconstruction. Even though the complete derivation of the turbulent GKS 
stress tensor is beyond the scope of this study, the idea is put forward that the turbulent 
GKS "activates" the high-order, non-linear terms as a response to insufficient numerical 
resolution and weighs them as a function of local rarefaction. In so doing, the turbulent 
GKS mimics the dynamics of more sophisticated turbulence models, such as the algebraic 
stress models. 

A link may be drawn to Chen's papei'^ which demonstrates that a turbulent gas-kinetic 
scheme may generate a turbulent stress tensor, mathematically consistent with non-linear 
turbulence models, that is, where the non-linear part can be expressed as follows: 
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r"^ = n^/k (c-^^^ + ^2 r + _^ c-i ——^ (47) 

* \ dxkdxk \dxkdxj dxk dxi J dxi dxj J ' 

Eq. |47]does not apply to the stress tensor generated by the GKS developed in this paper, 
as the latter is based on a first-order accurate Chap man- Enskog expansion (Eq. [9]); the 
fluxes' second-order accuracy is due to the analytical solution of the BGK model (Eq. |6]). 
Chen has proceeded from an expansion truncated to the second term and has not considered 
a discontinuous reconstruction. 

The turbulent GKS does not make any assumptions as to the type and behaviour of 
turbulence other than eddy dynamics can be represented as a "distribution of eddies". 
This makes it particularly elegant; nonetheless it also defines the limitations in practical 
cases. As a matter of fact, the turbulent GKS handles shock layers "naturally" - where 
rarefaction effects are significant, however, it cannot compensate for the lack of insight 
of two-equation models in critical situations such as in separated flow or in supersonic / 
hypersonic conditions. 



V. NUMERICAL IMPLEMENTATION 



The GKS fluxes indicated in Eqs. 26, 27 and 28 have been implemented in a finite 



volume steady-state solvei^nHSI^ relying on acceleration techniques such as 4-level multigrid 
and LU-SGS preconditioning, both implemented according to well-known paperi^^*^. The 
reconstruction include second and third order TVD/MUSCL schemes and fifth order WENO. 

The GKS has been implemented in both the directional splitting and multi-dimensional 
versions. In many flow cases, in accordance with May's finding^, no significant differences 
have been found. However, the multi-dimensional version provides better results in the 
presence of strong gradients, not aligned against grid directions. The stability of the multi- 
dimensional GKS also appears to be sensitive to various parameters such as CFL, limiters 
and cell size. Limiting primitive variables has proven more effective than limiting conser- 
vative variables, although its application cannot as yet be extended to the supersonic ramp 
flow. No attempts have been make to use multi-dimensional limiting. 

For comparison, results obtained from a Navier-Stokes solver are shown for most flow 
cases: this solver uses the same 4-level multigrid and LU-SGS preconditioning, whereas the 
advective fluxes are obtained from Roe's approximate Riemann solver and viscous fluxes 
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from central differences. 

The GKS usually remains stable at slightly higher CFL than the NSS, on average between 
4 and 10. However, the overall computational cost required by the GKS is higher: on 
average the GKS requires approximately twice the NSS time. The solver is sequential and 
no attempts have been made to optimize the code for speed. Turbulence has been modelled 
in both schemes with a standard linear k-u model, advanced in time in a segregated way. 
Artificial dissipation has been added in all GKS simulations. The value of C in Eqs. [19] and 



20 has been fixed heuristically, ranging from 0.25 to 0.75. The role of artificial dissipation 



is not yet completely clear; a small value is normally necessary to guarantee stability in 
the presence of a shock, whereas higher values may smooth the solution excessively or even 
affect the stability negatively. 

All the flow cases simulated concern transonic or supersonic flow at a Reynolds number 
ranging from one to several million. The computational meshes are stretched in order to 
place the first computational cell off solid walls within the laminar sublayer. 



VI. NUMERICAL EXPERIMENTS 

A. Transonic flow around a RAE 2822 airfoil in supercritical conditions 

The popular flow cases 9 and 10 investigated by Cook^^ around the RAE 2822 airfoil 
have been simulated with both the GKS and NSS. In case 9 the flow at M = 0.73 and 3.19° 
angle of attack, flow generates a strong normal shock at around 55 % of the airfoil chord 
without causing the separation of the boundary layer. In case 10, the flow at M = 0.745 
and 3.19° angle of attack, generates a stronger shock which causes an incipient separation, a 
more significant thickening of the boundary layer and a displacement upstream of the shock. 
The latter feature appears to be very challenging to capture for a solver based on a NSS 
and modelling turbulence with a standard, linear two-equation model. Reynolds number 
is slightly above 6 million in both cases. In the experiments, the boundary layer has been 
tripped at 3% of the airfoil chord on the upper and lower sides; the calculation is fully 
turbulent downstream of this point. Reconstruction is second-order TVD/MUSCL (with 
minmod limiters). 

Various C-type computational meshes have been used with a number of cells ranging from 
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625 X 125 to 928 x 160, with a cell-clustering at the expected shock position. The distance 
of the closest points from the airfoil surface is < 1, in wall units. Results obtained with 
the finest mesh are shown. The angles of attack have been reduced to take into account 
the wind tunnel corrections, to a = 2.80°. No special freestream boundary conditions have 
been applied to take into account the generation of vorticity. The distribution of pressure 
and skin friction coefficient obtained with the GKS and NSS for case 9 and 10 are shown in 
Fig. [T] and [2| respectively. In case 9 both schemes provide results which are in agreement 
with one another and with the experimental data. In case 10 the results obtained with the 
two schemes are very similar to one another, except for the shock area. The velocity profiles 
shown in Fig. [3] are reasonably close to each other and to experimental data. The results 
obtained from the GKS simulation are more in line with the esperimental data, in terms 
of position and thickness of the shock, whereas the NSS fails to capture the position of the 
shock correctly. The relatively poor prediction of separated fiows by the k-oj model in NS 
schemes is well-known and stems from the overprediction of turbulent shear stress in areas 
of incipient boundary layer separation: the errors found in this study from the NS scheme 
are in line with literature (refer for instance tc(^ and references therein). The prediction 
obtained from the GKS scheme shows an accuracy comparable to the one provided by more 
sophisticated, higher-order turbulence models in NSS. (Refer for instance to Wallin and 
Johansson'^ . 

Fig. |4] shows the relatively high values of the local turbulent Knudsen number (up to 
about 0.1) in the shock area and downstream, which is also the areas where the predictions 
provided by the two schemes differ most. 

B. Transonic flow around a NACA 0012 airfoil in supercritical conditions 

The NACA 0012 airfoil has been the object of several experimental investigations, includ- 
ing transonic and separated flow conditions. Two flow cases have been simulated with the 
multi-dimensional GKS and with the NSS: in the flrst one the flow at Mach M = 0.799 and 
angle of attack a = 2.86° causes an incipient separation of the boundary layer immediately 
downstream of the shock. In the second flow case, the flow at Mach M = 0.74 and angle 
of attack a = 4.86° causes the boundary layer to separate downstream of the shock and 
generate a large separated region. Reynolds number is 9 million in both cases. Both solvers 
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FIG. 1. Airfoil RAE2822, Case 9 Re = 6.3 x 10^, M = 0.73, angle of attack a = 3.19°. Pressure 
coefficient and skin friction. Continuous line: GKS, Dash-dotted line: NS, points: experimental 
data from CoolP^. 
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FIG. 2. Airfoil RAE2822, Case 10 Re = 6.2 x 10^, M = 0.745, angle of attack a = 3.19°. Pressure 
coefficient and skin friction. Continuous line: GKS, Dash-dotted line: NS, points: experimental 
data from CoolPl. 

have used the same C-type grid (for both cases) with 447 x 160 cells with 256 cells along 
the airfoil surface, with no cell clustering around the expected shock wave position. The 
distance of the closest points to the airfoil surface is < 1, in wall units. The angles of 
attack have been corrected to take into account the wind tunnel corrections, to a = 2.26° 
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FIG. 4. Airfoil RAE2822, Case 10. Contour of local turbulent Knudsen number on the airfoil 
upper side. Ten Mach lines have been added for reference. 

and a = 4.00° respectively. No special freestream boundary conditions have been applied 
to take into account the generation of vorticity. In the experiments, the boundary layer has 
been tripped at 5% of the airfoil chord on the upper and lower sides; the calculation is fully 
turbulent downstream of this point. Reconstruction is second-order TVD/MUSCL (with 
minmod limiters). 

Fig. [5] shows the pressure coefficient measured on the airfoil surface and calculated. 
In both cases, the GKS predictions indicate a shock position more in line with Harris' 
measurementpSI, as the Navier-Stokes solver with the linear k-uj model delays separation 
and fails to position the shock accurately in both cases, about 0.03 chords and 0.06 chords 
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respectively. 

However, both schemes fail to reproduce the correct behaviour of the pressure coefficient 
in the separated areas. This inaccurate prediction is due to the unsuitability of two-equation 
turbulence models to account for the behaviour of separated flow. Fig. [6] shows the local 
turbulent Knudsen number, calculated according to Eq. [STJ Inside the shocklayer, where 
GKS and Navier-Stokes solutions differ most, the turbulent Knudsen number reaches values 
of about 0.1. 
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X X 

FIG. 5. NACA 0012 Airfoil, on the LHS: Re = 9.0 x 10^, M = 0.799, angle of attack a = 2.86°. On 
the RHS: Re = 9.0 x 10*^, M = 0.74, angle of attack a = 4.86°. Pressure coefficient. Continuous 
line: GKS, Dash-dotted line: NS, points: experimental data from Harrid^ol^ 




FIG. 6. NACA 0012 Airfoil, Re = 9.0 x 10^, M = 0.74, angle of attack a = 4.86°. Contour of 
local turbulent Knudsen number on the airfoil upper side. Twenty Mach lines have been added for 
reference. 
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C. Transonic flow in Delery bump channel 



Figures [7] to 10 show the resuhs obtained from the simulation of the transonic bump 



flow (Case C), experimentally investigated by Delerj^. This is another well-known bench- 
mark case, suitable for validating turbulence models, as flow case features a strong shock 
- boundary layer interaction leading to a large flow separation immediately downstream of 
the shock. The Mach 0.615 duct flow impacts a one-side mounted ramp-semicircular bump, 
reaching approximately Mach 1.45 before the normal shock. The shock - boundary layer 
interaction generates the typical lambda structure, with the separation starting at the foot 
of the flrst leg. It is well known that predicting the position of the separation point and the 
size of the separated area are challenging for standard, linear two-equation models as they 
tend to delay separation and underpredict the extension of the separation. It is also well 
known that more sophisticated turbulence models including higher-order, non-linear terms 
provide more accurate predictions. Refer for instance to Goldber^^H 

The simulations have also been conducted with the solver based on the Navier-Stokes 
equations. Both solvers have used the same structured grid with 584 x 384 cells. The 
distance of the closest points to the wall surface is < 1, in wall units. In order to match 
the experimental position of the shock, the outlet pressure is heuristically adjusted. The 
two solvers therefore use slightly different values of outlet pressure. The calculation is fully 
turbulent. Reconstruction is second-order TVD/MUSCL (with minmod limiters). 

Fig. [T] (on the left hand side) shows good agreement between the distribution of the 
pressure obtained from both solvers and the pressure measured. However, the skin friction, 
shown in Fig. [7] on the right hand side, indicates that the separation predicted by the GKS 
is much more signiflcant than the one predicted by the NSS. The position of the separation 
with respect to the shock can be better observed from Fig. |8] where a few streamlines are 
shown: the separation predicted by the GKS starts correctly at the foot of the flrst shock leg, 
whereas the NSS predicts a signiflcantly delayed separation. Both solvers fail in accurately 
predicting the position of the re-attachment point and the shape of the separated region: in 
both cases the bubble's shape is too elongated. Whereas the GKS predicts a re-attachment 
point too far downstream, the NSS prediction seems to be more accurate simply because 
the separation is delayed. This inaccurate behaviour can be explained by the unsuitability 
of the two-equation models to predict correct values for the turbulence quantities in areas 
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of separated flow. Fig. |9] compares the distribution of Mach number obtained from the two 
solvers. 



Fig. 10 shows the distribution of turbulent Knudsen number, which assumes high values 
(around 0.2) immediately downstream of the boundary layer separation. This is also the 
area where GKS prediction differs most from the one obtained from the NSS. 
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FIG. 8. Delery bump channel flow. Streamlines highlighting the separation. Forty pressure contour 
lines added. GKS simulation on the left hand side, NSS simulation on the right hand side. 
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FIG. 10. Delery bump channel flow. Two representation of local turbulent Knudsen number. On 
left hand side, a pressure pseudo-color diagram is over-imposed. 

D. Shock-separated supersonic turbulent boundary layer at a compression 
corner 

The supersonic flow on a 24 degrees compression ramp has been simulated with the 
directional-splitting version of the GKS.The investigations carried out by Settle^^^*^ in the 
seventies have shown a large boundary layer separation. Schemes derived from the Navier- 
Stokes equations closed by linear two-equation turbulence models usually underestimate the 
size of the separation and fail to move the shock wave downstream (refer to Goldber^^. 
As a consequence of the separation area being too small, the predicted pressure rise is less 
gradual than the one measured. Non-linear models may correct this behaviour, as shown 
for instance by Goldber^^. However, the position of the re-attachment point and in general 
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the turbulence level downstream the shock is predicted in a rather inaccurate way also by 
non-linear models. This behaviour was also observed by Menteil^. 

The simulations have been executed on a structured grid with 288 x 256 cells. The distance 
of the closest points from the wall surface is y"'" < 1, in wall units. The prediction obtained 
with the GKS are accurate with regard to the position of the shock wave and the separation 
point; NSS with linear two-equation turbulence models tend to place the shock too little 
upstream of the corner thus further delaying separation (refer to Goldber^^. However the 
prediction obtained from the GKS are inaccurate for at least two reasons; (i) the pressure 
increase is much more abrupt than in Settles' measurement and (ii) the reattachment point 
is positioned too far downstream and the separation area has a shape which is much more 
elongated than the one observed by Settles. This non-physical behaviour can be explained 
by the unsuitability of two-equation models to predict accurate levels of the turbulence 
quantities in separated flow. 

The GKS provides a sufficiently clean picture of the shock system, shown in Fig. 12 
with the aid of pressure contours and selected Mach- and streamlines. The local turbulent 



Knudsen number is shown in Figs. 13 and 14 very high values are observed inside the shock 
value, precisely where the behaviour of the GKS differs most from the NSS. 



o 

a. 




0.001 
0.0008 
0.0006 
0.0004 
0.0002 


-0.0002 
-0.0004 
-0.0006 



o o ' 






o 


o 






o / - 


o 




o 







o / 


o 


o / 


o 


o / 




o / 


c 


cP / 







2 

x/do 



FIG. 11. Compression corner M = 2.85. Pressure and skin friction coefficient. Continuous fine: 
GKS, Dash-dotted fine: NS, points: experimental data from Settled^. 
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FIG. 12. Compression corner M = 2.85. Shock system represented with forty pressure contour 
hnes. The shp hne is also clearly visible. 




FIG. 13. Compression corner M = 2.85. Local Knudsen number. Twenty pressure contour lines 
have been added. 

VII. CONCLUSIONS 

GKS provide more accurate prediction than NSS in both the laminar and turbulent 
regimes. Despite the fact that the implementation of GKS into efficient and stable solvers is 
still challenging, it has the potential for improvement, as acceleration techniques and variable 
reconstruction can be tailored to this type of scheme. GKS are developed to compensate 
insufficient grid resolution with additional physics. In the turbulent regime, GKS exploit 
the underlying gas-kinetic theory to generate a non-linear turbulent stress tensor, whose 
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FIG. 14. Compression corner M = 2.85. Streamlines highlighting the separation. Forty pressure 
contour lines have been added. A pressure pseudo-color diagram is over-imposed. 

hydrodynamic limit bears some similarities to algebraic stress models. The higher-order 
contributions to the numerical fluxes depend on the ability of grid and solver to resolve 
the flow gradients; they are activated by the discontinuities in left and right reconstructed 
values. Moreover, the turbulent GKS exploits a valuable piece of information: the effects of 
collisions (or eddies interactions) on the advective fluxes. 

The numerical experiments shown in this paper - flow cases involving significant shock - 
boundary layer interaction - show that the turbulent GKS can improve on the NSS predic- 
tions in challenging flows, at the outer boundary of the applicability region of two-equation 
turbulence models. The analysis of the numerical simulations provides an additional link 
to more sophisticated turbulence models; the known weakness of linear k-uj models, to de- 
lay post-shock separation, is corrected by the GKS, in a manner close to the corrections 
provided by algebraic stress or non-linear turbulence models. 

This paper does not propose the turbulent GKS as the next turbulence model, but em- 
phasizes the fact that the mechanisms operated by GKS can be exploited to generate more 
powerful turbulence models. As a matter of fact, the GKS approach is potentially more 
interesting for Large Eddy Simulation. In this approach to turbulence modelling, only "sub- 
grid" turbulence is modelled, i.e. the flow dynamics that is not resolved by the grid. One 
of the biggest challenges - and the main reason for the success of dynamic model^ - is to 
design a subgrid model able to adjust itself to different flow conditions and grid sizes. 
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